week 7

In [1]:
import warnings
warnings.filterwarnings('ignore')
In [2]:
! jupyter nbextension enable --py --sys-prefix widgetsnbextension
Enabling notebook extension jupyter-js-widgets/extension...
      - Validating: ok
In [3]:
import numpy as np
import pandas as pd
import datetime as dt
import scipy.stats as stats 
In [4]:
# import holoviews as hv
# hv.extension('bokeh', logo=False)
# from holoviews.operation.datashader import aggregate
# import geoviews as gv
# import cartopy.crs as ccrs
In [5]:
from bokeh.plotting import figure, output_notebook, show
from bokeh.tile_providers import STAMEN_TERRAIN
output_notebook()

import datashader as ds
import datashader.transfer_functions as tf
from datashader.colors import Hot, viridis
from datashader.utils import export_image
from datashader.bokeh_ext import InteractiveImage

from holoviews.operation.datashader import aggregate

from functools import partial
from IPython.core.display import HTML, display

%matplotlib inline
import matplotlib.pyplot as plt
Loading BokehJS ...
In [6]:
def process_file(file_name):
    
    print('process of file: %s' %file_name)
    
    # читаем исходный файл с данными
    df = pd.read_csv(file_name, sep=',')
    
    # в некоторые годы с 2014 столбцы назывались иначе и был пробел перед названием колонки, унифицируем
    df.columns = [x.strip() for x in df.columns]
    if 'vendor_id' in df.columns:
        df.rename(columns={'vendor_id':'VendorID'}, inplace=True)
    if 'pickup_datetime' in df.columns:
        df.rename(columns={'pickup_datetime':'tpep_pickup_datetime'}, inplace=True)
    if 'dropoff_datetime' in df.columns:
        df.rename(columns={'dropoff_datetime':'tpep_dropoff_datetime'}, inplace=True)
    if 'rate_code' in df.columns:
        df.rename(columns={'rate_code':'RatecodeID'}, inplace=True)
    if 'RateCodeID' in df.columns:
        df.rename(columns={'RateCodeID':'RatecodeID'}, inplace=True)
    
    # удаляем записи где длительность поездки равна 0
    df = df[df.tpep_pickup_datetime != df.tpep_dropoff_datetime]

    # удаляем записи где количество пассажиров равно 0
    df = df[df.passenger_count != 0]

    # удаляем записи где дальность поездки равна 0
    df = df[df.trip_distance != 0.0]

    # удаляем записи, не входящих в прямоугольник Нью-Йорка
    df = df[(df.pickup_longitude <= -73.70001) & (df.pickup_longitude >= -74.25559) ] 
    df = df[(df.pickup_latitude >= 40.49612) & (df.pickup_latitude <= 40.91553) ]
    
    # переводим время посадки пассажиров из формата object в datetime
    df['tpep_pickup_datetime'] = pd.to_datetime(df['tpep_pickup_datetime'])
    
    # переводим время посадки пассажиров из формата object в datetime
    df['tpep_dropoff_datetime'] = pd.to_datetime(df['tpep_dropoff_datetime'])
    
    # получаем время поездки для записи 
    df['trip_time'] = (df['tpep_dropoff_datetime'] - df['tpep_pickup_datetime']).dt.total_seconds()
    
    # удаляем минуты или секунды из времени посадки пассажиров
    df['tpep_pickup_datetime'] = df.tpep_pickup_datetime.values.astype('<M8[h]')
    
    df.drop(['VendorID', 'passenger_count', 'trip_distance', 'RatecodeID', 'store_and_fwd_flag', 'fare_amount',
             'extra', 'mta_tax', 'tip_amount', 'tolls_amount', 'improvement_surcharge', 'total_amount',
             'trip_time'], axis=1, inplace=True)
    

    # берем из файла с идентификаторми зон данные по 2500 зонам
    regions_df = pd.read_csv('regions.csv', sep=';')

    # добавляем в таблицу данные по идентификаторам зон посадки
    binsx = regions_df.west[::50].values
    binsy = regions_df.south[:50].values
    res_pickup = stats.binned_statistic_2d(df.pickup_longitude, df.pickup_latitude, values=df.tpep_pickup_datetime, 
                                    statistic='count', bins=[binsx, binsy], expand_binnumbers=True)
    df['region_pickup'] = (res_pickup.binnumber[0] - 1) * 50 + res_pickup.binnumber[1]
    
    # добавляем в таблицу данные по идентификаторам зон высадки
#     binsx = regions_df.west[::50].values
#     binsy = regions_df.south[:50].values
#     res_dropoff = stats.binned_statistic_2d(df.dropoff_longitude, df.dropoff_latitude, values=df.tpep_pickup_datetime, 
#                                     statistic='count', bins=[binsx, binsy], expand_binnumbers=True)
#     df['region_dropoff'] = (res_dropoff.binnumber[0] - 1) * 50 + res_dropoff.binnumber[1]
    
    # получаем дополнительный признак - зона города (из 265), файлы с описанием должны быть в папке ./taxi_zones
    # !!!(этот участок обрабывается по 15 мин из-за spatial join 
#     zones_df = gpd.read_file('./taxi_zones/taxi_zones.shp').to_crs({'init': 'epsg:4326'})
#     zones_df.drop(['OBJECTID', 'Shape_Leng', 'Shape_Area', 'zone', 'borough'], axis=1, inplace=True)
#     shply_points = [Point(x, y) for x, y in zip(df['pickup_longitude'], df['pickup_latitude'])]
#     res = gpd.sjoin(gpd.GeoDataFrame(crs={'init':'epsg:4326'}, geometry=shply_points), zones_df, how='left', op='within')
#     df['zone'] = res['LocationID'].values
    
    # специальный ряд для таблиц где нужно посчитать общее количество
    df['num'] = np.ones(len(df))
    df = df.set_index('tpep_pickup_datetime', drop=False)
    
    return df
In [7]:
# %%time
# df_orig = pd.read_csv('yellow_tripdata_2016-06.csv')
In [8]:
# df_main_to_save = df_main[['pickup_longitude', 'pickup_latitude', 'region_pickup']]
In [9]:
# df_main_to_save.to_csv('yellow_tripdata_2016-06_edited.csv')
In [10]:
%%time
df_main = pd.read_csv('yellow_tripdata_2016-06_edited.csv', parse_dates=['tpep_pickup_datetime'])
df_main.set_index('tpep_pickup_datetime', inplace=True)
Wall time: 38.4 s
In [12]:
# загружаем предсказания количества поездок в июне 2016
predictions = pd.read_csv('predictions_2016-06.csv')

loc = predictions['id'].str.split('_')

regions = []; dates = []; hours = []; forwards = []
for x in loc:
    regions.append(x[0])
    dates.append(x[1])
    hours.append(x[2])
    forwards.append(int(x[3]))
    
df = pd.DataFrame()
df['regions'] = regions
df['datetime']=  [date+' '+hour for date, hour in zip(dates, hours)]
df['predicts'] = forwards
df['datetime'] = [pd.datetime.strptime(row, '%Y-%m-%d %H') for row in df['datetime']]
df['y'] = predictions['y']
df.set_index('datetime', inplace=True)
In [13]:
def get_day_predicts(df, forward, region=None, date=None, hour=None):
    if region is not None:
        df = df[df['regions'] == str(region)]
    df = df[df['predicts'] == int(forward)]
    df.is_copy = False
    df.index = df.index + dt.timedelta(hours=1)
    if date is not None:
        df = df[date]
    if hour is not None:
        df = df[df.index.hour == hour]['y'].values
    return df
In [14]:
def get_bootstrap_samples(df_loc, predictions):
    regions = df.regions.unique().astype(int)
    result = {}
    for n_samples, region in zip(predictions, regions):
        n_samples = int(round(n_samples, 0))
        data = df_loc[df_loc['region_pickup'] == region]
        if n_samples != 0: 
            index = np.random.choice(data.shape[0], size=int(n_samples))
            row = data.iloc[index, [0,1]].values
            result[region] = row
        else:
            result[region] = [0, 0]
    pred_coords = pd.DataFrame()
    for x in result.values():
        y = pd.DataFrame(x)
        pred_coords = pd.concat([pred_coords, y])
#     print(pred_coords.columns)
    pred_coords.columns = ['pickup_longitude', 'pickup_latitude']
    pred_coords.index = list(range(len(pred_coords)))
    return pred_coords
In [15]:
def get_real_df(date, hour, region=None):
    df_real = df_main[['pickup_longitude', 'pickup_latitude', 'region_pickup']]
    if date is not None:
        df_real = df_real[date]
    if hour is not None:
        df_real = df_real[df_real.index.hour == hour]
    df_real = df_real[df_real['region_pickup'].isin(df.regions.unique().astype(int)) ]
    if region is not None:
        df_real = df_real[df_real['region_pickup'] == region]
    df_real.drop('region_pickup', axis=1, inplace=True)
    df_real.index = list(range(len(df_real)))
    return df_real
In [16]:
def latlng_to_meters(df, lat_name, lng_name):
    lat = df[lat_name]
    lng = df[lng_name]
    origin_shift = 2 * np.pi * 6378137 / 2.0
    mx = lng * origin_shift / 180.0
    my = np.log(np.tan((90 + lat) * np.pi / 360.0)) / (np.pi / 180.0)
    my = my * origin_shift / 180.0
    df.loc[:, 'pickup_x'] = mx
    df.loc[:, 'pickup_y'] = my
In [17]:
import param, paramnb
In [18]:
df_loc = pd.read_csv('time_region_102_pickup_count.csv')
df_loc = df_loc.set_index('datetime')
df_loc = df_loc['2016-06-01':]
In [19]:
class GetData(param.Parameterized):
    hour = param.Integer(default=12, bounds=(0,24))
    date = param.ObjectSelector(default='2016-06-01', 
       objects=np.unique([x.strftime('%Y-%m-%d') for x in df_loc.index.astype('<M8[D]')]), 
                check_on_set=True)

1) карты с визуализацией реального и прогнозируемого спроса на такси в выбираемый пользователем момент времени

In [20]:
paramnb.Widgets(GetData, button=True, next_n=4)
In [21]:
date = GetData.date
hour = GetData.hour
# получаем координаты реальных проездок в выбранный день и час
df_real = get_real_df(date, hour, region=None)
# получаем предсказанное количество поездок в выбранный день и час (из ранее полученной модели) 
if (date == '2016-06-30') & (hour > 18):
    prediction_nums = get_day_predicts(df=df, forward=1, region=None, date=date, hour=hour)
    print(prediction_nums.shape)
    prediction_nums2 = get_day_predicts(df=df, forward=6, region=None, date=date, hour=hour)
    print(prediction_nums2.shape)
    prediction_nums = pd.concat([pd.DataFrame(prediction_nums), pd.DataFrame(prediction_nums2[-5:])], axis=0)
    print(prediction_nums.shape)
    # переиндексируем таблицу
    hours = pd.date_range('2016-06-01 00:00:00', periods=720, freq='H')
    prediction_nums.index = hours
else:
    prediction_nums = get_day_predicts(df=df, forward=1, region=None, date=date, hour=hour)
# бутстрапируем коордианты из реальной выборки поячеечно для предсказанного количеста поездок в каждой ячейке
df_pred = get_bootstrap_samples(df_main, prediction_nums)
# добавляем в файл данные в метрах для получения карт ниже
latlng_to_meters(df_real, 'pickup_latitude', 'pickup_longitude')
latlng_to_meters(df_pred, 'pickup_latitude', 'pickup_longitude')
In [31]:
# NY = x_range, y_range = ((df_real['pickup_x'].min(), df_real['pickup_x'].max()), 
#                           (df_real['pickup_y'].min(), df_real['pickup_y'].max()))
NYC = x_range, y_range = ((-8242000,-8210000), (4965000,4990000))

plot_width  = 900
plot_height = 900

def base_plot(tools='pan, box_zoom, wheel_zoom, save, reset', 
              plot_width=plot_width, plot_height=plot_height, **plot_args):
    p = figure(tools=tools, plot_width=plot_width, plot_height=plot_height,
        x_range=x_range, y_range=y_range, outline_line_color=None,
        min_border=0, min_border_left=0, min_border_right=0,
        min_border_top=0, min_border_bottom=0, **plot_args)
    
    p.axis.visible = False
    p.xgrid.grid_line_color = None
    p.ygrid.grid_line_color = None
    return p

def create_image_pred(x_range, y_range, w=plot_width, h=plot_height):
    cvs = ds.Canvas(plot_width=w, plot_height=h, x_range=x_range, y_range=y_range)
    agg = cvs.points(df_pred, 'pickup_x', 'pickup_y')
    img = tf.shade(agg.where(agg>np.percentile(agg,90)), cmap=viridis, how='eq_hist')
    return tf.dynspread(img, threshold=0.5, max_px=4)

def create_image_real(x_range, y_range, w=plot_width, h=plot_height):
    cvs = ds.Canvas(plot_width=w, plot_height=h, x_range=x_range, y_range=y_range)
    agg = cvs.points(df_real, 'pickup_x', 'pickup_y')
    img = tf.shade(agg.where(agg>np.percentile(agg,90)), cmap=viridis, how='eq_hist')
    return tf.dynspread(img, threshold=0.5, max_px=4)
In [32]:
# карта с визуализацией реального спроса
export = partial(export_image, export_path="export")
p = base_plot()
p.add_tile(STAMEN_TERRAIN)
export(create_image_real(*NY), "NY1")
InteractiveImage(p, create_image_real)
Out[32]:
In [33]:
# карта с визуализацией предсказанного спроса
export = partial(export_image, export_path="export")
p = base_plot()
p.add_tile(STAMEN_TERRAIN)
export(create_image_pred(*NY), "NY2")
InteractiveImage(p, create_image_pred)
Out[33]:

2) временной ряд фактического и прогнозируемого спроса на такси в выбираемой области.

In [25]:
class GetCell(param.Parameterized):
    cell = param.ObjectSelector(default=df.regions.unique().astype(int)[0], 
                                objects=df.regions.unique().astype(int), check_on_set=True)
In [26]:
from sklearn.metrics import mean_absolute_error
In [27]:
paramnb.Widgets(GetCell, button=True, next_n=4)
In [34]:
region = GetCell.cell
# получаем предсказанное количество поездок в выбранный день и час (из ранее полученной модели) 
prediction_nums = get_day_predicts(df=df, forward=1, region=region, date=None, hour=None)
# т.к. в файле предсказания на 6 часов вперед, то последние 5 записей не вошли в вышеукзанный файл, добавляем их 
prediction_nums2 = get_day_predicts(df=df, forward=6, region=region, date=None, hour=None)
prediction_nums = pd.concat([prediction_nums, prediction_nums2[-5:]])
# переиндексируем таблицу
hours = pd.date_range('2016-06-01 00:00:00', periods=720, freq='H')
prediction_nums.index = hours
pred_arr = prediction_nums['y'].values
In [35]:
df_loc = pd.read_csv('time_region_102_pickup_count.csv')
df_loc = df_loc.set_index('datetime')
df_loc = df_loc['2016-06-01':]
real_arr = df_loc.loc[:,str(region)].values
In [36]:
plt.figure(figsize=(18,6))
plt.plot(pred_arr, 'r', label='predicted')
plt.plot(real_arr, 'b', label='real values')
plt.title('2016 June, cell: {}, MAE: {}, XGBRegressor'.format(region, mean_absolute_error(real_arr, pred_arr)))
plt.ylabel('pickups')
plt.xlabel('hours in june 2016')
plt.legend()
plt.show()
In [ ]: